#include "HY1C_out.h"
/* ============================================================================ */
/* module l1_modis_hdf.c - functions to read MODIS L1B for MSL12                */
/*                                                                              */
/* Written By: B. Franz, NASA/SIMBIOS, January 2003.                            */
/*                                                                              */
/* ============================================================================ */

#include "l1_modis_hdf.h"
#include "l12_proto.h"
#include "hdf_utils.h"

#define EOSMETALEN 32768

#define NEVBANDS 15
#define NEMBANDS 16
#define MAXBANDS 16
#define NDET     10

static int have_rescaled = 0;
static int have_extract  = 0;

static int32 sd_id;
static int32 sd_id_g;
static int32 spix   = 0;
static int32 sscan  = 0;  /* start scan in original granule */


/* RSMAS corrections to Terra IR radiances (20,22,23,27,28,29,31,32) */
static float radoff[NBANDSIR][NDET] = {
{-0.000972, 0.014200, 0.000022, 0.000238,-0.000003, 0.003760, 0.002337, 0.000084, 0.008848, 0.005050},
{-0.002300, 0.000600, 0.000236,-0.000280,-0.000003, 0.002798, 0.003496, 0.018035, 0.002942, 0.001787},
{-0.003833, 0.003657, 0.002118, 0.000798,-0.000003, 0.000208, 0.000399, 0.000553, 0.000258, 0.021128},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{-0.000423,-0.000242,-0.000330,-0.000065,-0.000001,-0.000006, 0.000064, 0.000398, 0.000362, 0.000322},
{-0.000433,-0.000246,-0.000222,-0.000148,-0.000001, 0.000068, 0.000221, 0.000221, 0.000437, 0.000303}
};
static float mfact[NBANDSIR] = {0.0,0.0,0.0005,0.0,0.0,0.0,0.0028,0.00278};

/* RSMAS AQUA modis corrections for some bands and detectors */
static float ch20d2cor  =  0.996666;    /* channel 20 detector 2  */
static float ch20d3cor  =  1.0044444;   /* channel 20 detector 3  */
static float ch20d10cor =  1.0033333;   /* channel 20 detector 10 */
static float ch22d1cor  =  1.0030075;   /* channel 22 detector 1  */
static float ch22d8cor  =  1.0007518;   /* channel 22 detector 8  */
static float ch22d9cor  =  1.0015037;   /* channel 22 detector 9  */
static float ch22d10cor =  1.0022556;   /* channel 22 detector 10 */
static float ch23d1cor  =  1.0192857;   /* channel 23 detector 1  */


/* ----------------------------------------------------------------------------------- */
/* compute_alpha() - computes the frame rotation angle for polarization correction     */
/*                                                                                     */
/* algorithm provided by: F. S. Patt, SAIC, February 2003.                             */
/* ----------------------------------------------------------------------------------- */
void compute_alpha(float lon[], float lat[], float senz[], float sena[], double mnorm[3],
                   int32 npix, float alpha[])
{
    static double radeg = 57.29577951;

    int32 ip, i;

    double slon, clon;
    double slat, clat;
    double szen, czen;
    double sazi, cazi;

    double e[3], n[3], v[3], r[3], s[3];
    double sdvcr, vdr, sdr, sdv;

    /* invert mirror normal */
    for (i=0; i<=2; i++)
      r[i] = -mnorm[i];

    for (ip=0; ip<npix; ip++) {

      slon = sin(lon [ip]/radeg);
      clon = cos(lon [ip]/radeg);
      slat = sin(lat [ip]/radeg);
      clat = cos(lat [ip]/radeg);
      szen = sin(senz[ip]/radeg);
      czen = cos(senz[ip]/radeg);
      sazi = sin(sena[ip]/radeg);
      cazi = cos(sena[ip]/radeg);

      /* pixel coordinate system (north, east, vertical) in ECR */

      e[0] = -slon;
      e[1] =  clon;
      e[2] =  0.0;

      n[0] = -slat * clon;
      n[1] = -slat * slon;
      n[2] =  clat;

      v[0] =  clat * clon;
      v[1] =  clat * slon;
      v[2] =  slat;

      /* sensor view vector in ECR */

      for (i=0; i<=2; i++)
          s[i] = e[i] * szen * sazi
               + n[i] * szen * cazi
               + v[i] * czen;

      /* compute rotation angle (alpha) from pixel normal (v) to mirror */
      /* normal (r) about sensor view vector (s)  (Wertz, p. 757)       */

      sdvcr = s[0] * (v[1]*r[2]-v[2]*r[1])
            + s[1] * (v[2]*r[0]-v[0]*r[2])
            + s[2] * (v[0]*r[1]-v[1]*r[0]);

      vdr   = v[0]*r[0] + v[1]*r[1] + v[2]*r[2];

      sdr   = s[0]*r[0] + s[1]*r[1] + s[2]*r[2];

      sdv   = v[0]*s[0] + v[1]*s[1] + v[2]*s[2];

      /* negated to be consistent with Gordon et al. */

      alpha[ip] = -(radeg * atan2(sdvcr,(vdr-sdr*sdv)) - 90.0);

    }

}


/* ----------------------------------------------------------------------------------- */
/* get_hdfeos_meta() - retrieves meta-data strings from HDF EOS core or archive meta.  */
/*                                                                                     */
/* B. Franz, SAIC, March 2003.                                                         */
/* G. Fireman, SAIC, March 2010: revised                                               */
/* ----------------------------------------------------------------------------------- */
int32 get_hdfeos_meta(int32 sd_id, char *metastr, char *name, char *data) {
  char  eosmeta[EOSMETALEN] = "";
  char *p = NULL;
  char *q = NULL;
  char *data2;

  data2 = strdup(data);

  if (SDreadattr(sd_id, SDfindattr(sd_id,metastr), (VOIDP) eosmeta) != 0) {
    HY1C_out("Error reading meta string %s\n",metastr);
    return(-1);
  }
  data [0] = '\0';
  data2[0] = '\0';

  /* find object */
  if ((p = strstr(eosmeta,name)) == NULL)
    return(-1);

  /* extract value */
  if ((p = strstr(p,"VALUE")) == NULL)
    return(-1);
  if ((p = strstr(p,"=")) == NULL)
    return(-1);
  p++;
  if ((q = strstr(p,"\n")) == NULL)
    return(-1);
  strncpy(data,p,q-p);
  data[q-p] = '\0';

  /* remove enclosing quotes as needed */
  if (strstr(data,"\"")){
    data2 = strdup(data);
    p = strstr(data2,"\"");
    p++;
    q = p + strlen(data2);
    if (strstr(p,"\""))
      q = strstr(p,"\"");
    strncpy(data,p,q-p);
    data[q-p] = '\0';
    free(data2);
  }

  return(0);
}


/* ----------------------------------------------------------------------------------- */
/* G. Fireman, SAIC, March 2010 */
/* ----------------------------------------------------------------------------------- */
int32 get_modis_orbit(int32 sd_id, filehandle *file) {

    char result[80] = "";
    int  yy,mm,dd,hh,mn,ss, year,day,msec;
    char frac[10] = "";

    yy = mm = dd = hh = mn = ss = 0;

    get_hdfeos_meta(sd_id,"CoreMetadata.0","ORBITNUMBER",result);
    file->orbit_number = (int32_t) atoi(result);

    get_hdfeos_meta(sd_id,"CoreMetadata.0","EQUATORCROSSINGLONGITUDE",result);
    file->orbit_node_lon = atof(result);

    get_hdfeos_meta(sd_id,"CoreMetadata.0","EQUATORCROSSINGDATE",result);
    sscanf(result, "%d-%d-%d", &yy, &mm, &dd);

    get_hdfeos_meta(sd_id,"CoreMetadata.0","EQUATORCROSSINGTIME",result);
    sscanf(result, "%d:%d:%d.%s", &hh, &mn, &ss, frac);

    ymdhms2ydmsec(yy,mm,dd,0,0,0, (int32_t*) &year, (int32_t*) &day, (int32_t*) &msec);
    sprintf(file->node_crossing_time,"%04d%03d%02d%02d%02d%s",year,day,hh,mn,ss,frac);

    return(0);
}


/* ----------------------------------------------------------------------------------- */
/* B. Franz, SAIC, December 2007.                                                      */
/* ----------------------------------------------------------------------------------- */
int32 get_modis_calfile(int32 sd_id,char *file) {

   char  eosmeta[EOSMETALEN] = "";
    char *p = NULL;
    char *q = NULL;
    int   i;

    if (SDreadattr(sd_id, SDfindattr(sd_id,"CoreMetadata.0"), (VOIDP) eosmeta) != 0) {
        HY1C_out("Error reading meta string CoreMetadata.0\n");
        return(-1);
    }

    if ((p = strstr(eosmeta,"INPUTGRANULE")) == NULL)
      return(-1);
    if ((p = strstr(p,"VALUE")) == NULL)
      return(-1);
    if ((p = strstr(p,"=")) == NULL)
      return(-1);
    p += 3;
    if ((q = strstr(p,")")) == NULL)
      return(-1);
    q -= 1;
    strncpy(file,p,q-p);
    file[q-p] = '\0';

    p = file;
    for (i=0; i<strlen(file); i++) {
        if (file[i] != '\\' && file[i] != '\"' && file[i] != 'S') {
          if ((q = strstr(&file[i],"UR_")) != &file[i])
            *p++ = file[i];
          else
            i+= 2;
        }
    }
    *p = '\0';

    return(0);
}


/* ----------------------------------------------------------------------------------- */
/* openl1_modis_hdf() - opens a MODIS L1B file for reading.                            */
/*                                                                                     */
/* B. Franz, SAIC, February 2003.                                                      */
/* ----------------------------------------------------------------------------------- */
int openl1_modis_hdf(filehandle *file)
{
    int32  npix;
    int32  nscan;

    int32  sds_id;
    int32  rank;
    int32  dims[3];
    int32  type;
    int32  numattr;

    int32  itemp;
    uint32 gflags[8];


    /* Open the HDF input file */
    sd_id = SDstart(file->name, DFACC_RDONLY);
    if(sd_id == FAIL){
        fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
        __FILE__,__LINE__,file->name,DFACC_RDONLY);
        return(HDF_FUNCTION_ERROR);
    }

    /* Get pixel and scan dimensions */
    sds_id = SDselect(sd_id,SDnametoindex(sd_id,("EV_1KM_RefSB")));
    if (SDgetinfo(sds_id,NULL,&rank,dims,&type,&numattr) == -1) {
        fHY1C_out(stderr,"-E- %s line %d: error getting dimension info.\n",
        __FILE__,__LINE__);
        return(HDF_FUNCTION_ERROR);
    }
    npix  = dims[2];
    nscan = dims[1];

    /* Do we have rescaled band data to extend the R, G, B and NIR bands? */
    if (getHDFattr(sd_id,"Rescaled Ocean R","EV_250_Aggr1km_RefSB",(VOIDP)&itemp) == 0) {
        if(want_verbose)
            HY1C_out("File contains rescaled RGBN data.\n");
        have_rescaled = 1;
    }

   /* Is this an L1B made from an L1A extract */
    if (getHDFattr(sd_id,"Extract Pixel Count","",(VOIDP)&itemp) == 0 && itemp>0) {
        /* The L1B and GEO is scan extracted, but padded to full scan */
        npix  = itemp;
        getHDFattr(sd_id,"Extract Pixel Offset","",(VOIDP)&itemp);
        spix  = itemp;
        getHDFattr(sd_id,"Extract Line Offset","",(VOIDP)&itemp);
        sscan = itemp*NDET;
        have_extract = 1;
        HY1C_out("File was generated from L1A extract of size %d x %d\n",(int)npix,(int)nscan);
        HY1C_out("  Pixels %d - %d\n",(int)(spix+1), (int)( spix+ npix));
        HY1C_out("  Lines  %d - %d\n",(int)(sscan+1),(int)(sscan+nscan));
    }

    file->npix   = npix;
    file->ndets  = NDET;
    file->nscan  = nscan;
    file->sd_id  = sd_id;
    get_modis_calfile(sd_id,file->calfile); /* read list of inputs to L1B */


    /* Open the HDF geolocation input file.*/
    sd_id_g = SDstart(file->geofile, DFACC_RDONLY);
    if(sd_id_g == FAIL){
        HY1C_out("Error opening geolocation file.\n");
        fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
            __FILE__,__LINE__,file->geofile,DFACC_RDONLY);
        return(HDF_FUNCTION_ERROR);
    }

    /* Has geolocation been corrected for terrain height? */
    if (getHDFattr(sd_id_g,"Cumulated gflags","",(VOIDP)&gflags) == 0) {
      file->terrain_corrected = (int32_t)(gflags[5]==0);
    } else HY1C_out("-E- %s line %d: Error reading gflags.\n",__FILE__,__LINE__);

    get_modis_orbit(sd_id_g, file); /* read orbit info */


    return(LIFE_IS_GOOD);
}


/* ----------------------------------------------------------------------------------- */
/* readl1_modis_hdf() - reads 1 line (scan) from a MODIS L1B file, loads l1rec.        */
/*                                                                                     */
/* B. Franz, SAIC, February 2003.                                                      */
/* ----------------------------------------------------------------------------------- */
int readl1_modis_hdf(filehandle *file, int32 scan, l1str *l1rec)
{
    static int   firstCall = 1;
    static float pi = 3.141592654;

    static uint16 *data = NULL;

    static uint16 *data_250 = NULL;
    static uint16 *data_500 = NULL;
    static  int16 *hgt  = NULL;

    //needed
    static  int16 *senz = NULL;
    static  int16 *sena = NULL;
    static  int16 *solz = NULL;
    static  int16 *sola = NULL;

    /* Mapping from EV_1KM block to bands of interest   */
    /* MODIS band number = 8,9,10,11,12,13lo,14lo,15,16 */
    /* Wavelength = 412,443,488,531,551,667,678,748,870 */

    static int   bmap[NBANDS]   = {0,1,2,3,4,5,7,9,10};
    static float m   [NEVBANDS];
    static float b   [NEVBANDS];

    static float refl_250_m[2];
    static float refl_250_b[2];
    static float refl_500_m[5];
    static float refl_500_b[5];

    static float Fo_blue=205.733;
    static float Fo_green=183.76;
    static float Fo_red=157.807;
    static float Fo_ir=97.174;

    /* Mapping from EM_1KM block to bands of interest */
    /* MODIS band number = 20,22,23,27,28,29,31,32    */
    /* Wavelength = 3.7, 3.9, 4, 6.7, 7.3, 8.5, 11,12 */
    /* 20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36*/

    static int   bmapir[NBANDSIR] = {0,2,3,6,7,8,10,11};
    static float m_ir  [NEMBANDS];
    static float b_ir  [NEMBANDS];

    static float m_cirrus;
    static float b_cirrus;

    static float64 lastTAIsec = -1;
    static float64 TAIsec, usec, dsec;
    //needed
    static int16  year, day;
    static double fsol;
    static int16   mside, badmside;
    static double mnorm[3];
    static int32  lastframe = -1;

    int32 npix   = (int32)file->npix;
    int32 nbands = (int32)file->nbands;
    int32 nbandsir = NBANDSIR;
    int32 frame  = scan / NDET;
    int32 detnum = scan % NDET;
    int32 iw,ip,ipb;

    static float *Fobar;
    double esdist;

    if (firstCall) {

        firstCall = 0;

        /* We need this for the reflectance to radiance conversion */
        rdsensorinfo(file->sensorID,file->input->evalmask,"Fobar", (void **) &Fobar);

        if ((data = (uint16 *) calloc(npix*MAXBANDS, sizeof(uint16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating data space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"reflectance_scales","EV_1KM_RefSB",(VOIDP)m) != 0) {
             HY1C_out("-E- %s line %d: Error reading reflectance scale attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"reflectance_offsets","EV_1KM_RefSB",(VOIDP)b) != 0) {
            HY1C_out("-E- %s line %d: Error reading reflectance offset attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"radiance_scales","EV_1KM_Emissive",(VOIDP)m_ir) != 0) {
            HY1C_out("-E- %s line %d: Error reading radiance scale attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"radiance_offsets","EV_1KM_Emissive",(VOIDP)b_ir) != 0) {
            HY1C_out("-E- %s line %d: Error reading radiance offset attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((data_250 = (uint16 *) calloc(npix*2, sizeof(uint16))) == NULL) {
          HY1C_out("-E- %s line %d: Error allocating data_250 space.\n",
                 __FILE__,__LINE__);
          return(1);
        }

        if ((data_500 = (uint16 *) calloc(npix*2, sizeof(uint16))) == NULL) {
          HY1C_out("-E- %s line %d: Error allocating data_500 space.\n",
                 __FILE__,__LINE__);
          return(1);
        }

        if (getHDFattr(sd_id,"reflectance_scales","EV_250_Aggr1km_RefSB",(VOIDP)refl_250_m) != 0) {
          HY1C_out("-E- %s line %d: Error reading reflectance scale attribute.\n",
                 __FILE__,__LINE__);
          return(1);
        }

        if (getHDFattr(sd_id,"reflectance_offsets","EV_250_Aggr1km_RefSB",(VOIDP)refl_250_b) != 0) {
          HY1C_out("-E- %s line %d: Error reading reflectance offset attribute.\n",
                 __FILE__,__LINE__);
          return(1);
        }

        if (getHDFattr(sd_id,"reflectance_scales","EV_500_Aggr1km_RefSB",(VOIDP)refl_500_m) != 0) {
          HY1C_out("-E- %s line %d: Error reading reflectance scale attribute.\n",
                 __FILE__,__LINE__);
          return(1);
        }

        if (getHDFattr(sd_id,"reflectance_offsets","EV_500_Aggr1km_RefSB",(VOIDP)refl_500_b) != 0) {
          HY1C_out("-E- %s line %d: Error reading reflectance offset attribute.\n",
                 __FILE__,__LINE__);
          return(1);
        }


        if (getHDFattr(sd_id,"reflectance_scales","EV_Band26",(VOIDP)&m_cirrus) != 0) {
             HY1C_out("-E- %s line %d: Error reading reflectance scale attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if (getHDFattr(sd_id,"reflectance_offsets","EV_Band26",(VOIDP)&b_cirrus) != 0) {
             HY1C_out("-E- %s line %d: Error reading reflectance offset attribute.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((hgt = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating hgt space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((solz = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating solz space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((sola = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating sola space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((senz = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating senz space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

        if ((sena = (int16 *) calloc(npix, sizeof(int16))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating sena space.\n",
                   __FILE__,__LINE__);
            return(1);
        }

    }

    if (frame != lastframe) {

        /* Get mirror side (0=A,1=B) and mirror normal (in ECR) */
        READ_SDS_ID(sd_id_g,"Mirror side",&mside, frame,0,0,0, 1,1,1,1);
        READ_SDS_ID(sd_id_g,"T_inst2ECR" ,&mnorm, frame,0,0,0, 1,3,1,1);

        /* Fix mirror-side if set to -1 in geolocation (used as array index) */
        if (mside < 0) {
            mside = 0;
            badmside = 1;
        } else
            badmside = 0;

        /* Read granule time and compute scan times */
        READ_SDS_ID(sd_id_g,"EV start time",(float64 *)&TAIsec, frame,0,0,0, 1,1,1,1);

        /* Work-around for single-scan time errors */
        if (TAIsec < 0.0) {
            HY1C_out("-W- %s: bad time in geolocation file at frame %d, using previous.\n",
                   __FILE__,(int)frame);
            TAIsec = lastTAIsec;
        } else
            lastTAIsec = TAIsec;

        usec = TAIsec + 725846400.0;
        unix2yds(usec,&year,&day,&dsec);

        *(l1rec->year) = (int32_t)year;
        *(l1rec->day)  = (int32_t)day;
        *(l1rec->msec) = (int32_t)(dsec*1000.0);

        /* get earth-sun distance correction for this frame */
        esdist = esdist_(l1rec->year,l1rec->day,l1rec->msec);
        fsol   = pow(1.0/esdist,2);

        lastframe = frame;
    }

    *(l1rec->year) = (int32_t)year;
    *(l1rec->day)  = (int32_t)day;
    *(l1rec->msec) = (int32_t)(dsec*1000.0);

    /* Get position and path geometry */
    READ_SDS_ID(sd_id_g,"Longitude",l1rec->lon,    scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"Latitude" ,l1rec->lat,    scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"Height"         ,hgt,     scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"SolarZenith"    ,solz,    scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"SolarAzimuth"   ,sola,    scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"SensorZenith"   ,senz,    scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"SensorAzimuth"  ,sena,    scan,spix,0,0, 1,npix,1,1);

    for (ip=0; ip<npix; ip++) {
        l1rec->height[ip] = (float)hgt[ip];
        l1rec->solz[ip]   = (float)solz[ip] * 0.01;
        l1rec->sola[ip]   = (float)sola[ip] * 0.01;
        l1rec->senz[ip]   = (float)senz[ip] * 0.01;
        l1rec->sena[ip]   = (float)sena[ip] * 0.01;

        l1rec->pixnum[ip] = spix + ip;

        if (l1rec->lon[ip] < -181.0 || l1rec->lon[ip] > 181.0 ||
            l1rec->lat[ip] <  -91.0 || l1rec->lat[ip] >  91.0 ||
            badmside)
          l1rec->navfail[ip] = 1;
    }

    /* Compute polarization frame rotation angles */
    compute_alpha(l1rec->lon, l1rec->lat, l1rec->senz, l1rec->sena,
                  mnorm, npix, l1rec->alpha);
    if (mside == 0) {
        if ((file->input->evalmask & MSKMODMIR1) != 0)
            for (ip=0; ip<npix; ip++) l1rec->navfail[ip] = 1;
    } else {
        if ((file->input->evalmask & MSKMODMIR2) != 0)
            for (ip=0; ip<npix; ip++) l1rec->navfail[ip] = 1;
    }

    /* Read L1B VNIR data, scale to radiance, and copy relevant bands to L1 record  */
    /* MODIS format is band sequential.  We want band interlaced.  We use imbedded  */
    /* scale factors to get reflectance, then convert to radiance through mean F0.  */
    /* Note that use of MCST radiance scale factor gives different result.          */

    /* If data is available, saturated bands 10,12,13,and 16 will be extended using */
    /* a rescaled version of the data from the 250 and 500-meter bands              */
    /*  2 (EV_1KM_RefSB)                0 (EV_500_Aggr1km_RefSB)  Blue              */
    /*  4 (EV_1KM_RefSB)                1 (EV_500_Aggr1km_RefSB)  Green             */
    /*  5 (EV_1KM_RefSB)                0 (EV_250_Aggr1km_RefSB)  Red               */
    /* 10 (EV_1KM_RefSB)                1 (EV_250_Aggr1km_RefSB)  NIR               */

    READ_SDS_ID(sd_id,"EV_1KM_RefSB",data, 0,scan,spix,0, NEVBANDS,1,npix,1);

    READ_SDS_ID(sd_id,"EV_250_Aggr1km_RefSB",data_250, 0,scan,spix,0, 2,1,npix,1);
    READ_SDS_ID(sd_id,"EV_500_Aggr1km_RefSB",data_500, 0,scan,spix,0, 2,1,npix,1);

    /*
    for (ip=0; ip<npix; ip++)
      HY1C_out("%d %d\n",mside,l1rec->navfail[ip]);
    */

    for (iw=0; iw<nbands; iw++) {

        l1rec->Fo[iw] = Fobar[iw] * fsol;

        for (ip=0;ip<npix; ip++) {

            ipb = ip*NBANDS+iw;
            l1rec->Lt[ipb] = BAD_FLT;

            /* skip night data for visible bands */
            if (l1rec->solz[ip] >= SOLZNIGHT)
                continue;

            /* check for sentinel values and flag as appropriate */
              if (data[bmap[iw]*npix+ip] >= 65500) {
                  switch (data[bmap[iw]*npix+ip]) {
                case 65527: /* not in earth view */
                     l1rec->navfail[ip] = 1;
                    break;
                case 65521: /* UNRESCALED_HIGH_SI */
                case 65520: /* RESCALED_SI */
                  l1rec->hilt[ip] = 1;
                  l1rec->Lt[ipb] = 1000.0;
                  if (have_rescaled) {
                    switch (bmap[iw]) {
                    case 2: /* Blue */
                      l1rec->Lt[ipb] = (data_500[0*npix+ip] - refl_500_b[0])
                        * refl_500_m[0] * l1rec->Fo[iw]/pi;
                      break;
                    case 4: /* Green */
                      l1rec->Lt[ipb] = (data_500[1*npix+ip] - refl_500_b[1])
                        * refl_500_m[1] * l1rec->Fo[iw]/pi;
                      break;
                    case 5: /* Red */
                      l1rec->Lt[ipb] = (data_250[0*npix+ip] - refl_250_b[0])
                        * refl_250_m[0] * l1rec->Fo[iw]/pi;
                      break;
                    case 10: /* NIR */
                      l1rec->Lt[ipb] = (data_250[1*npix+ip] - refl_250_b[1])
                        * refl_250_m[1] * l1rec->Fo[iw]/pi;
                      break;
                    }
                  } else {
                    switch (bmap[iw]) {
                    case 2: /* Blue */
                      l1rec->Lt[ipb] = (data_500[0*npix+ip] - refl_500_b[0])
                        * refl_500_m[0] * Fo_blue/pi;
                      break;
                    case 4: /* Green */
                      l1rec->Lt[ipb] = (data_500[1*npix+ip] - refl_500_b[1])
                        * refl_500_m[1] * Fo_green/pi;
                      break;
                    case 5: /* Red */
                      l1rec->Lt[ipb] = (data_250[0*npix+ip] - refl_250_b[0])
                        * refl_250_m[0] * Fo_red/pi;
                      break;
                    case 10: /* NIR */
                      l1rec->Lt[ipb] = (data_250[1*npix+ip] - refl_250_b[1])
                        * refl_250_m[1] * Fo_ir/pi;
                      break;
                    }
                    /* Fix for old files with no aggr info */
                    if (l1rec->Lt[ipb] <= 0.0)
                       l1rec->Lt[ipb] = 1000.0;
                  }
                  break;
                case 65533: /* saturation */
		case 65529:
                  l1rec->hilt[ip] = 1;
                  l1rec->Lt[ipb] = 1000.0;
                  if (have_rescaled == 0) {
                    switch (bmap[iw]) {
                    case 2: /* Blue */
                      l1rec->Lt[ipb] = (data_500[0*npix+ip] - refl_500_b[0])
                        * refl_500_m[0] * Fo_blue/pi;
                      break;
                    case 4: /* Green */
                      l1rec->Lt[ipb] = (data_500[1*npix+ip] - refl_500_b[1])
                        * refl_500_m[1] * Fo_green/pi;
                      break;
                    case 5: /* Red */
                      l1rec->Lt[ipb] = (data_250[0*npix+ip] - refl_250_b[0])
                        * refl_250_m[0] * Fo_red/pi;
                      break;
                    case 10: /* NIR */
                      l1rec->Lt[ipb] = (data_250[1*npix+ip] - refl_250_b[1])
                        * refl_250_m[1] * Fo_ir/pi;
                      break;
                    }
                    /* Fix for old files with no aggr info */
                    if (l1rec->Lt[ipb] <= 0.0)
                       l1rec->Lt[ipb] = 1000.0;
                    break;
                  }
                default:    /* other stuff */
                  break;
                }

            } else
              l1rec->Lt[ipb] = (data[bmap[iw]*npix+ip] - b[bmap[iw]])
                * m[bmap[iw]] * l1rec->Fo[iw]/pi;
          }
    }

    /* Read L1B IR data, scale to radiance, and copy relevant bands to L1 record */

    READ_SDS_ID(sd_id,"EV_1KM_Emissive",data, 0,scan,spix,0, NEMBANDS,1,npix,1);

    for (iw=0; iw<nbandsir; iw++) {

        for (ip=0;ip<npix; ip++) {

            ipb = ip*NBANDSIR+iw;
            l1rec->Ltir[ipb] = 0.0;

            /* check for sentinel values and flag as appropriate */
              if (data[bmapir[iw]*npix+ip] >= 65500) {
              ;
            } else {
                l1rec->Ltir[ipb] =
                    (data[bmapir[iw]*npix+ip] - b_ir[bmapir[iw]]) * m_ir[bmapir[iw]]/10.;

                /* apply RSMAS radiance corrections */
                if (l1rec->sensorID == MODIST) {
                    l1rec->Ltir[ipb] *= (1.0-radoff[iw][detnum]);
                    if (mside == 1) l1rec->Ltir[ipb] *= (1.0 + mfact[iw]);
                }
                if (((file->input->evalmask & SSTMODS) != 0) &&
                    (l1rec->sensorID == MODISA || l1rec->sensorID == HMODISA)) {
                    if (iw == 0 && detnum == 2-1) {
                        l1rec->Ltir[ipb] *= ch20d2cor;
                    } else if (iw == 0 && detnum == 3-1) {
                        l1rec->Ltir[ipb] *= ch20d3cor;
                    } else if (iw == 0 && detnum == 10-1) {
                        l1rec->Ltir[ipb] *= ch20d10cor;
                    } else if (iw == 1 && detnum == 1-1) {
                        l1rec->Ltir[ipb] *= ch22d1cor;
                    } else if (iw == 1 && detnum == 8-1) {
                        l1rec->Ltir[ipb] *= ch22d8cor;
                    } else if (iw == 1 && detnum == 9-1) {
                        l1rec->Ltir[ipb] *= ch22d9cor;
                    } else if (iw == 1 && detnum == 10-1) {
                        l1rec->Ltir[ipb] *= ch22d10cor;
                    } else if (iw == 2 && detnum == 1-1) {
                        l1rec->Ltir[ipb] *= ch23d1cor;
                    }
                }
            }
          }
    }

    /* Read cirrus band */
    READ_SDS_ID(sd_id,"EV_Band26",data, scan,spix,0,0, 1,npix,1,1);

    for (ip=0;ip<npix; ip++) {
        l1rec->rho_cirrus[ip] = 0.0;
          if (data[ip] < 65500) {
          l1rec->rho_cirrus[ip] = (data[ip] - b_cirrus) * m_cirrus / cos(l1rec->solz[ip]/RADEG);
        }
    }

    l1rec->sensorID = file->sensorID;
    l1rec->npix     = file->npix;
    l1rec->detnum   = (int32_t)detnum;
    l1rec->mside    = (int32_t)mside;

    /* Convert IR bands to brightness temperature */
    radiance2bt(l1rec,-1);

    return(LIFE_IS_GOOD);
}


/* ----------------------------------------------------------------------------- */
/* readl1_lonlat_modis_hdf() - reads 1 line (scan) from a MODIS L1B file, loads  */
/*                             loads only the lon and lat values l1rec.          */
/*                                                                               */
/* Don Shea, SAIC, May 2009.                                                     */
/* ----------------------------------------------------------------------------- */
int readl1_lonlat_modis_hdf(filehandle *file, int32 scan, l1str *l1rec)
{
    int32 npix   = (int32)file->npix;

    /* Get position and path geometry */
    READ_SDS_ID(sd_id_g,"Longitude",l1rec->lon, scan,spix,0,0, 1,npix,1,1);
    READ_SDS_ID(sd_id_g,"Latitude" ,l1rec->lat, scan,spix,0,0, 1,npix,1,1);

    return(LIFE_IS_GOOD);
}


int closel1_modis_hdf(filehandle *file)
{
   if ((sd_id_g != sd_id) && SDend(sd_id_g)) {
      fHY1C_out(stderr,"-E- %s line %d: SDend(%d) failed for file, %s.\n",
      __FILE__,__LINE__,sd_id,file->geofile);
      return(HDF_FUNCTION_ERROR);
   }
   if (SDend(sd_id)) {
      fHY1C_out(stderr,"-E- %s line %d: SDend(%d) failed for file, %s.\n",
      __FILE__,__LINE__,sd_id,file->name);
      return(HDF_FUNCTION_ERROR);
   }

   return(LIFE_IS_GOOD);
}



/* L1B DN* Sentinel Values:       */
/* L1A_SCAN_DATA_MISSING_SI 65535 */
/* L1A_DN_MISSING_SI        65534 */
/* SATURATED_DETECTOR_SI    65533 */
/* ZERO_POINT_DN_SI         65532 */
/* DEAD_DETECTOR_SI         65531 */
/* RSB_DN_STAR_BELOW_MIN_SI 65530 */
/* TEB_OR_RSB_GT_MAX_SI     65529 */
/* AGGREGATION_FAIL_SI      65528 */
/* SECTOR_ROTATION_SI       65527 */
/* TEB_B1_NOT_CALCULATED    65526 */
/* NAD_CLOSED_UPPER_SI      65500 */

/* Bit QA Flags int32 1=True; 0=False
bit  0 Moon in SV Port
bit  1 Spacecraft Maneuver
bit  2 Sector Rotation
bit  3 Negative Radiance Beyond Noise Level
bit  4 PC Ecal on
bit  5 PV Ecal on
bit  6 SD Door Open
bit  7 SD Screen Down
bit  8 SRCA On
bit  9 SDSM On
bit 10 Outgassing
bit 11 Instrument Standby Mode
bit 12 Linear Emissive Calibration
bit 13 DC Restore Change
bit 14 BB/Cavity Temperature Differential
bit 15 BB Heater On
bit 16 Missing Previous Granule
bit 17 Missing Subsequent Granule
bit 18 SRCA calibration mode
bit 19 SRCA calibration mode
bit 20 Moon within the SV keep-out box for Reflective Solar bands
bit 21 Moon within the SV keep-out box for Emissive bands
bit 22 All space-view data bad for any Reflective Solar band
bit 23 All blackbody data bad for any Reflective Solar band
bit 24 Dropped scan(s) between previous granule and granule being calibrated.
bit 25 Dropped scan(s) between  granule being calibrated and subsequent granule.
bit 26
*/
